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Abstract 

Linear kinetic Monte Carlo particle transport models are frequently 
employed in fusion plasma simulations to quantify atomic and surface 
effects on the main plasma flow dynamics. Separate codes are used for 
transport of neutral particles (incl. radiation) and charged particles (trace 
impurity ions). Integration of both modules into main plasma fluid solvers 
provides then self consistent solutions, in principle. The required inter- 
faces are far from trivial, because rapid atomic processes in particular 
in the edge region of fusion plasmas require either smoothing and re- 
sampling, or frequent transfer of particles from one into the other Monte 
Carlo code. We propose a different scheme here, in which despite the 
inherently different mathematical form of kinetic equations for ions and 
neutrals (e.g. Fokker-Planck vs. Boltzmann collision integrals) both types 
of particle orbits can be integrated into one single code. We show that 
the approximations and shortcomings of this "single sourcing" concept 
(e.g., restriction to explicit ion drift orbit integration) can be fully tol- 
erable in a wide range of typical fusion edge plasma conditions, and be 
overcompensated by the code-system simplicity, as well as by inherently 
ensured consistency in geometry (one single numerical grid only) and (the 
common) atomic and surface process modules. 

This is a pre peer reviewed version which has been submitted to Com- 
puter Physics Communications 

1 Introduction 

A common computational approach for scrape off layer plasmas of fusion exper- 
iments arc combined packages which consist of two coupled codes: a fluid solver 
for the main plasma (deuterium and tritium) and a kinetic solver for neutrals 
(deuterium, tritium, tungsten, carbon...), see e.g. [2 [21 ISl In addition a 
kinetic code for impurity ions (tungsten, carbon, beryllium...) can be added, 
e.g. [5]. This approach avoids improper assumptions of the fluid modeling of 
impurity ions, which are instantaneous thermalization of the impurity particles 
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and the lack of kinetic effects as a whole. Many standalone kinetic Monte Carlo 
impurity transport codes have been developed in the past, e.g. [6l [71 HI |9] . In 
general coupling of codes requires averaging of involved exchanged quantities 
due to different coordinate systems and grids used by the codes resulting in 
additional inaccuracies which degrade the output. Thus direct integration and 
streamlining of codes is favorable. The present approach combines kinetic treat- 
ment of neutrals with the kinetic treatment of impurity ions in linear approx- 
imation, where the plasma background is fixed. The kinetic neutral transport 
code EIRENE is supplemented by a trace ion module (TIM), which comprises 
the numerical methods for solving a linear drift kinetic equation with Monte 
Carlo approach. The main difference of the TIM to other impurity codes is that 
it is fully part of EIRENE and many routines for describing the neutral parti- 
cle dynamics are reused for treating charged particles. Thus EIRENE is now 
capable of solving three types of equations in one single code: the Boltzmann 
Equation for the neutrals, [Ullll], a Boltzmann like equation describing photon 
gas transport problems, |13] . and TIM, which now allows solving the linear drift 
kinetic equation for impurity ions. 

This paper is organized as follows: The current status of the Monte Carlo 
transport code EIRENE is summarized in the next chapter. Then the basic 
idea of the trace ion module is given in chapter 3. The linear drift kinetic 
model is shortly summarized in chapter 4 and its numerical implementation in 
the EIRENE code is described in chapter 5. The scope of application of the 
present trajectory integration method is analyzed in chapter 6 and the collision 
operator is verified in chapter 7. Finally the extended code is applied to a 
realistic modeling scenario for the MAST tokamak, where basic physics features 
of the trace ion module are verified. In addition the simulation results have 
been compared with CCD camera images for two different wave lengths. The 
summary is given in chapter 7. Technical details arc described in the appendix. 



2 Overview of the Monte Carlo Code EIRENE 

A brief overview of the current status of the Monte Carlo transport code EIRENE 
pn [T2] is given. This code package has been developed for modeling neutral 
gas and radiation transport problems in magnetically confined plasmas. It is 
a multi-species code solving simultaneously a system of time dependent or sta- 
tionary linear kinetic transport equations of almost arbitrary complexity. Linear 
kinetic transport problems are characterized by a given background plasma dis- 
tribution function /^(v), which is usually reconstructed from fluid quantities. 
There are no restrictions concerning the geometric complexity of the problem, 
in general any discretization from 0-3 dimensional computational domains is 
possible. It has been mostly used together with fluid transport codes, where it 
supplies the necessary particle, momentum and energy sources from the neu- 
tral particle dynamics to the fluid codes. In return these fluid codes provide 
the plasma background for EIRENE. Surface processes (physical and chemical 
sputtering as well as reflection) arc treated by EIRENE with sputter models 
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Figure 1: Outline of the Monte Carlo code EIRENE and where the trace ion 
module is interfacing this code 



and reflection data bases, see [M]. Atomic and molecular data is obtained from 
external data bases, most notably ADAS [15]. EIRENE further has access to 
the HYDKIN cross section data base for hydrocarbon molecules, Refs. [23l [24]. 
for simulation of the complicated catabolism mechanisms of these hydrocarbons 
in the fusion plasma. Since EIRENE solves a general kinetic equation radiation 
transport problems (photon gas simulations) can be treated as well, one such 
example is modeling of High-Intensity Discharge (HID) lamps with EIRENE 
P"5] . The code structure of the EIRENE code and possibly other Monte Carlo 
codes for solving linear transport problems is outlined in figure [1] An approx- 
imate model for transport of ionized particles along magnetic field lines has 
previously also been part of the EIRENE code. This model is now expanded 
by the trace ion module, which treats kinetic transport physics of ions with less 
simplifications. 
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3 Basic Idea of the Trace Ion Module 



The trace ion module is an extension of the EIRENE code for solving a hnear 
drift kinetic equation, see e.g. |18| , for describing transport of charged particles 
in the existing code framework. The idea was to just add new orbit following 
routines for ions as well as a Fokker-Planck collision operator. In contrast to 
the neutral particle dynamics which is governed by the Boltzmann equation the 
drift kinetic equation is of Fokker-Planck type. Boltzmann equations describe 
discontinuous jump processes, where short range interactions between particles 
are dominant. On the other hand Fokker-Planck equations describe diffusion 
processes. The sample paths are described by stochastic ordinary differential 
equations, referred to as Langevin equations. Random sampling as well as time 
discretization is required for numerically solving stochastic ordinary differential 
equations. For more details on stochastic ordinary differential equations see e.g. 
[HI [19]. 

A combination of discontinuous jump processes and time discrete integra- 
tion of stochastic ordinary differential equations is possible whenever the time 
discrete integration algorithm can be formulated as a jump process. In this case 
the Langevin equations which correspond to the drift kinetic equation have been 
solved by splitting the particle motion into a deterministic part which is inter- 
rupted by an artificial collision event. This collision event takes into account 
the effective action of the Coulomb force on the particle trajectory during one 
time step as well as the change of the velocity vector due to the electromagnetic 
forces in guiding center approximation. The major problem is that EIRENE is 
working with distances rather than time steps for calculating trajectories of the 
neutrals. In fact EIRENE determines the flight distance of a neutral particle 
from an effective mean free path and samples the type of collision at the point 
of collision from a discrete distribution of probabilities of the involved collision 
processes. Only once this distance and the point of collision is known a time 
step can be calculated, which is simply the amount of time necessary to reach 
the point of collision at the current particle speed. 

The natural coordinate system for following of guiding center orbits is aligned 
with the magnetic field, which decouples parallel and perpendicular motion and 
significantly simplifies trajectory integration. This is conflicting with the neu- 
tral particle dynamics which is treated in Cartesian coordinates. In practice 
the EIRENE code is working with unstructured grids in a Cartesian coordi- 
nate system and the integration of the guiding center equations has to fit into 
this scheme. An obvious choice for orbit integration are Runge Kutta methods, 
which are frequently used in many kinetic plasma transport codes. In principle 
such methods could be incorporated into EIRENE under the expense of addi- 
tional particle tracking on the grid. If such methods are applicable in EIRENE 
without a significant performance loss is to be investigated. Computationally 
less demanding but highly accurate methods for orbit integration are Adams 
Bashforth backward methods. Currently a 4-step Adams-Bashforth backward 
formula has been incorporated into EIRENE, which makes the calculation of the 
effective guiding center velocity simple and efficient. Moreover numerical errors 
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due to Cartesian coordinates have been mitigated by significantly improving 
the interpolation and differentiation techniques on the numerical grid. Interpo- 
lation methods from finite element theory have been applied, where grid cells 
possess shape functions allowing to interpolate and differentiate electric and 
magnetic field vectors locally in each grid cell. Currently three types of cells are 
supported: 3-node triangles, 4-node quadrangles and 4-node tetrahedrons. 

Coulomb collisions are treated as field collisions, where test particles interac- 
tion with the background plasma is described by evaluating Trubnikov/Rosenbluth 
potentials This is done one the fly for non-Maxwellian distribution func- 
tions, which enables thermal force effects on a kinetic level. Since this collision 
operator in guiding center coordinates is singular an implicit method has been 
developed, which handles the singularity occurring in the perpendicular drift 
coefficient. The present approach, the details of which are summarized in |F1 
shows a solution for overcoming this problem in the framework of a Monte Carlo 
approach, while solutions for the same problem in context with finite difference 
methods have already been reported in e.g. [55]. The location of the new orbit 
integration and collision routines of the trace ion module in the code structure 
of EIRENE is sketched in figure [H 

4 Some remarks on the linear Drift 
Kinetic Model 

A summary of the drift kinetic model is given. Drift kinetic theory applies 
whenever the Larmor radius of the considered ion is much smaller than the 
gradient length of the external magnetic fields, ^ B/\WB\. The Larmor 
radius for ions with mass m and charge number Z is defined by = ^2/^^, 
where V2 is the perpendicular velocity component and the cyclotron frequency 
is given by O = ZeB/m. The magnetic field vector is denoted by B, its absolute 
value by B and the unit vector by b. The charge of the electron is e. The 
following notations are applied, 

vi =v - h, E'li = E • b, 

where the velocity component parallel to the magnetic field is denoted with Vi . 
The drift kinetic equation for stationary magnetic and electric fields and a fixed 
plasma background (linear approximation) can then be written as 

d ^ 

^ iv2f) = -Vy(y • V2f) - dk {[Vk + Ak]v2f) 

fe=i 
1 ^ 

+ 2 E ^fc^' iDkiV2f) + S, (1) 

which governs the time evolution of the distribution function f{t, y, vi, W2). For 
a detailed derivation of this equation see e.g. |18| . This equation is already 
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written in Fokker-Planck form which allows straight forward construction of 
appropriate Monte Carlo methods. Therein the guiding center velocity y and 
the time derivatives of the reduced phase space velocities wi, V2 are defined by 
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This approach covers the drifts due to radial electric field and magnetic field 
inhomogeneities as well as the mirror effect and the electrostatic force. Coulomb 
collisions are represented by the drift and diffusion coefhcients, A), and Dki 
respectively. A comprehensive report on evaluating these coefficients for non- 
Maxwellian plasmas can be found in jlOj . Therein appropriate drift and diffusion 
coefficients in the reduced (wi, V2) phase space are given, 

= ,A^, 

OVi 

av2 2 ^2 dv2 
D22 = A-^, 



The Trubnikov/Roscnbluth potentials (f> and ijj arc specialized for taking into 
account friction as well as the thermal force effect. The derivation of the explicit 
form of these potentials is given in |B1 which differs from the derivation given in 
[TU] in the following respects: a different method has been used for integration of 
the potential functions and the differentials required for calculating the thermal 
force effect have been expanded to cover the perpendicular direction, which will 
enable thermal force effects in radial direction in the future. It is noted that the 
notations and definitions of the potentials derived in llOj have been reused in 
this work. The standard Monte Carlo approach for diffusion in real space with 
a constant diffusion coefficient is used, where the diffusion matrix is of the form 

Dij = d^/dyidyj {eiekD±V2f) ■ (4) 
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Therein the vector e is defined to be perpendicular to the magnetic field in 
radial direction lying in the poloidal plane. The factor A is given by 

where and denotes particle density and charge of the background particles. 
The Coulomb logarithm A replaces a singular integral which results from the 
infinite range of the electrostatic potential. This integral is usually cut off at 
a certain impact parameter which is of the order of the Debye radius. The 
Coulomb logarithm has been taken constant through out the whole work. As 
a typical value for fusion relevant plasmas A = 13.5 has been used. Since 
impurities in the SOL plasma consist of many different species and each of these 
species can have more than one charge state additional sources/sinks represented 
by S in equation [1] account for ionization, recombination {Si, Sr) and external 
particle creation and annihilation processes (Q^)- "^^^ complete source/sink 
term reads 

S = -Sf{v2ff + Sf-\v2ff-' 



5 Numerical Implementation 

Random walks of neutral particles in the EIRENE code are constructed by 
evaluating the distance d to the next collision according to the effective local 
mean free path of the involved collision processes. Once this distance is known 
EIRENE updates the position of a particle according to 

Ynew = Void + evd (5) 

where particle positions are denoted by y and the unit velocity vector is e.u . At 
the new position the next collisional event is executed and the velocity vector 
V ■ ey of the particle is instantaneously changed. This process is continued until 
a time limit is reached or the particle is absorbed at a surface. Macroscopic 
quantities are estimated from the random walks of the particles on the compu- 
tational grid by applying a track length estimator. Random walks of this type, 
referred to as Markov chains, solve the Boltzmann equation for the neutrals in 
terms of distribution functions. 

In the present approach the dynamics of ionized particles is governed by 
stochastic ordinary differential equations corresponding to [TJ The first order 
Euler-Maruyama Method for discretization of these equations is given by 

Xt+At ^Xt+ AiXt)At + B(Xt)VAtC, (6) 

where X = {Xt,t > 0} is the state vector of the system, which is represented 
in this case by the position of the particle y and the reduced (ui,W2) velocity 
phase space. Local drift coefficients are denoted hy A= (y, wi.2+^1.2), diffusion 
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coefficients by B and C is a normal distributed random number with mean and 
variance 1. Diffusion coefficients B can be constructed from the corresponding 
diffusion matrix {D = BB^) in[TJ which in principal is a five dimensional matrix 
covering diffusion in real space as well as in velocity space. An appropriate 
method for constructing B in velocity space is given in [E] The initial state of a 
new born particle, due to e.g. a surface process, is the initial position and the 
parallel and perpendicular velocity components of the 3D velocity vector, which 
EIRENE assigns to newly created particles. 

For incorporating the discretized Langevin equations into the Monte Carlo 
framework of the EIRENE code the actual algorithm has to be formulated as 
a jump process. In particular the position of the ionized particle has to be 
updated in the same way as for the neutrals according to eq. (O, but with an 
effective guiding center velocity 

«^l|y!l e,„^y/||y||. 

This effective guiding center velocity vector is constructed at each point of 
collision and takes into account the effect of the electromagnetic force as well as 
the Coulomb force both acting on the particle over a time period At. Currently 
a 4-step Adams Bashforth backward formula is used for evaluating y, the details 
of which are given in [A] For calculating a distance to the next Coulomb collision 
event, 

d^v Atcc, (7) 

an appropriate time step Atcc has to be supplied. This is not necessarily 
the time step used for the actual collision process, but it allows EIRENE to 
determine the shortest distance for the next type of collision. Coulomb collisions 
have to be executed in any case, no matter what other collision event is executed. 
In practice the time step for Coulomb collisions and furthermore the distance to 
the next Coulomb collision has to be much smaller than the mean free path of 
the other involved collision processes for proper resolving the thermalization of 
the test particles with the plasma background. In fact Atcc is just the upper 
limit for the actual time step. The actual time step for updating the guiding 
center position can be calculated after the flight distance d has been determined 
by EIRENE from 

At - d/||y|| < Atcc- (8) 

For resolving thermalization processes Atcc can be adjusted as a fraction, e.g. 
1/100, of the local equilibration time. 

The actual algorithm works in the following manner: after the distance to 
the next collision has been determined by EIRENE and the time step At is 
known the particle position is updated according to equation [5j where e^, • d is 
replaced with At ■y{t). In the next step post collision phase space velocities are 
calculated by 

^r>o.t-col ^ ^re-col ^ ^ ^^ocH ^ (9) 

where normalized quantities are used, = ''Ji with the inverse thermal velocity 
a ~ ^j2Tb/mb. This is an effective collision where in addition to the Coulomb 
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interaction, Ax'**""^, also the acceleration due to the electric field and the mirror 
force, Ax'''^*, is accounted for. The latter are calculated at the current particle 
position y post-col s-nd are given by 

Axf* = a(\v^v.2^\ At. (10) 

It is noted that consistency requires to evaluate the deterministic velocity in- 
crements Ax''*^* in the same way as y. The stochastic velocity increments are 
evaluated in any case according to 

2 

^^StOch ^ ^^^^ ^ J2 B^.VMCj, (11) 

with two random numbers Ci.2- 

The singularity of this collision operator is expressed in the perpendicular 
drift coefficient A2, which is singular for X2 ~^ 0. If friction forces only are 
considered a different set of coordinates might have allowed to get around this 
difliculty, but for a numerical implementation of the thermal force effect in 
the kinetic Monte Carlo approach the present coordinates (wi,W2) are the most 
simple and computationally the most efficient ones. The present algorithm 
has been supplemented by an implicit part, which treats slow particles (with 
respect to 172) below a certain threshold. The additional computational expense 
is about 10%, where it is noted that without this implicit part the energy of 
the test particles is overestimated and conservation of energy of the collision 
operator is severely violated. Details are described in |f1 

The next step of this algorithm is finding interpolated values and derivatives 
of the magnetic and electric fields at the point of collision Ypost-coU 

■ ( post — col post — col -n / \ T71 / \\ 

y Xi ,X2 ^^{y post — col) T post — col/ j ■ 

and evaluating the new effective guiding center velocity by applying an appro- 
priate integration method. Finally the particle position is advanced to the next 
point of collision. In this approach xi £md X2 are actually dummy variables 
which are only required for constructing the guiding center velocity. In fact 
an ionized particle is described in EIRENE by the position y and the velocity 
y in the six dimensional 3D in real and 3D in (guiding center) velocity space. 
It is noted that in 2D simulations the toroidal contribution of the curvature 
drift cannot be resolved. The magnetic field is given on the poloidal plane only, 
which is in fact a cylindrical approximation. A possible axis symmetric toroidal 
magnetic field component has the form B(r, 0) = BqRq/ Re^p, where is the 
unit vector in toroidal direction. The toroidal part of the curvature drift is then 
recovered as 

VCD = -b X (b X (V X b)) = --ej,, (12) 

R 
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where is a unit vector pointing upward in the direction of the axis of the 
tokamak (the plane in which the torus is located is the e^: x e^-planc). The 
direction of this drift is aligned with the VB drift, which is assumed to point 
downwards here. Especially for spherical tokamaks, e.g. MAST, the toroidal 
curvature is one order of magnitude higher as for other tokamaks and the effect 
of the toroidal curvature has a big impact on modeling results. 

6 Guiding Center Orbit Integration 

The quality of orbit integration in the presented algorithm is determined by 
the evaluation of the effective guiding center velocity at each time step. If y is 
calculated from the local plasma and magnetic field parameters at the current 
simulation time the resulting jump process is identical to the first order Eu- 
ler method and fits perfectly to EIRENE by lowest computational costs. The 
drawback of applying the Euler method for guiding center orbit integration in 
fusion plasmas is the accumulating numerical error, which leads to an artificial 
outward drift of the ions. This is only tolerable for short living particles where 
closed orbits are not to be expected. For example the life time of C"*""*" in the 
divertor region of the MAST tokamak is 200/zs according to ionization rates 
of ADAS, [H]. This is long enough to fully thermalize. On the other hand the 
bounce time estimated for a circular magnetic field (major radius i?o chosen for 
MAST edge plasma) and a particle at a thermal speed of 10*m/s is ti, = 10~^s. 
The lifetime of C++ is two orders of magnitude lower than the time for passing 
at least a single orbit, which justifies to use the simple Euler orbit integrator in 
this case. 

For treating particle motion of particles with longer life times more accurate 
orbit integration methods have been sought. In fusion research Runge Kutta 
methods are frequently used for guiding center orbit following codes. One such 
example is ASCOT, [71 15]. Implementing &x\ n~ th order Runge Kutta meth- 
ods into EIRENE requires the evaluation of n intermediate steps, for which 
each time the whole geometry module (locate particle on grid, interpolate value 
inside grid cell) has to be called. Since the geometry module is the most de- 
manding part of EIRENE with respect to computing power such methods are 
not suitable. Promising methods, which allows to improve the orbit integration 
of EIRENE while keeping the additional computational effort low, are Adams 
Bashforth backward formulas. The effective velocity vector needed for advanc- 
ing the particle position is constructed from previous velocity values rather 
than from calculating intermediate steps, which requires almost no additional 
computing power just a negligible increase of storage. The details of the 4- 
step backward method which has been incorporated into the EIRENE code is 
described in The reason why a 4-step backward method has been used is 
because backward formulas taking two and three steps were observed to be un- 
stable. Nevertheless these formulas are required to start the particle orbit. This 
unstable behavior might be due to extreme time step differences, which can 
occur in EIRENE. In particular in front of surfaces the time step can change 
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e.g. from At = 10"'^ to At = 10" because EIRENE is moving a particle from 
one cell to the cell face of the next one for keeping track of the particles on the 
grid. This method can be easily extended to 5 or more step method, but this 
has not been considered yet. 

As an illustrative example of the performance of the backward orbit integra- 
tor is shown in figure [H The evolution of the radial coordinate at the mid 
plane positive, outer zero crossing of the orbit) of a closed collision less or- 
bit is plotted vs. time. The orbit, which has been calculated with the Euler 
method, continuously drifts outward while the radial coordinate of the orbit 
obtained with the 4 step Adams-Bashforth keeps constant (300 orbits transits 
are shown). The outward drift per orbit is not completely removed with the 
Adams-Bashforth method, but compared to the Euler integrator it is consid- 
erably lower. For this particular orbit, which has been integrated with a step 
size of At = 10~^s, the radial outward drift per orbit has been found to be 
dr Adams/ dr ^ lO^^m/r ^ drEuici/dr ^ lO^'^m/r, where r denotes the amount 
of time to pass one full orbit. 




Figure 2: Time evolution of radial coordinate at mid plane of closed banana orbit 
for simple first order orbit integration and 4-stcp Adams Bashforth backward 
method, magnetic field in circular approximation with dimensions of MAST 
closed field line region, collision less particle launched close to scparatrix 



7 Verification of the Coulomb Collision Opera- 
tor 

Whereas the steady state behavior of the present Cou- lomb collision operator 
has been well tested in Ref. [10], verification of the dynamics of the collision 
operator is missing. For this reason the time evolution of corresponding mo- 
ments of the kinetic distribution function have been estimated and the results 
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have been confronted to available analytical relaxation time approximations. 
In particular, slowing down of particles, temperature equilibration and temper- 
ature isotropization have been checked. The most simple estimator has been 
used for calculating moments of the distribution function, M, which is given 
M = l/(AFiiVi) ^ ■ g{vi), where g{vi) is any function of Vi, e.g. particle energy 
mvf/2. A slab case for the EIRENE code has been set up existing of one cell 
with periodic boundary conditions and uniform magnetic field in z direction. 
Typical plasma parameters for fusion edge plasmas have been adjusted in these 
simulations, ~ lOeV, nu+ = lO^^m"^ and u = Oms~^. A point source of 
test particles in the cell center has been placed, where ions with isotropic, 
monoenergetic velocity distribution have been released at an initial energy of 
leV. The time interval for the simulations was l/2ms, where Uniuai = /xs and 
t final = 500 /is and the time step has been fixed to At = 10~^s. An analytic 
expression for the equilibration of plasmas at different temperatures but no rel- 
ative drift velocity is given by dT/dt = {Tt — T) where the characteristic 
frequency is 



Temperature isotropization of a plasma with constant overall temperature but 
different Ty and Tj. is described by dAT/dt ^ -uiAT where AT = Ti - T2 and 
the characteristic frequency 



where 7 = 1/^2 + (3 - A) atanh {VA)/Va) and yl = 1 - T\\o/Tj_o. In this 



case particles have to be sampled from a bi-Maxwellian distribution function. 
The overall temperature, T = (Tj| -|- 2Tl)/3, is constant in this case and energy 
is solely transferred from parallel to perpendicular direction. Slowing down 
means that the macroscopic flow velocity of the test particles adjusts to the 
flow velocity of the plasma background, described by a shifted Maxwellian. The 
governing relaxation equation is d{v) = —vsi{v) — Ufa) where the slowing down 
frequency is 



The slowing down of test particles has been checked for a beam of particles 
aligned with the magnetic field. Results of the Monte Carlo simulations as 
well as the solutions of the ordinary differential equations for slowing down, 
thermalization and isotropization have been plotted in figure |3l Therein the 
time evolution of T and u, normalized to background plasma temperature as 
well as flow velocity, and the time evolution of AT is shown. Analytical results 
and simulation are in good agreement for all three processes. 

Moreover the present collision operator includes the thermal force effect, the 
numerical implementation of which has been adapted for the EIRENE code 




(13) 




(14) 



i^s = -Ma^0o(x)/x- 



(15) 
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Figure 3: Comparison of relaxation time approximations and Monte Carlo sim- 
ulation, left ordinate corresponds to normalized (v) ju (slowing down, blue line) 
and T/Tb (thermalization, black line) and the right ordinate corresponds to AT 
(isotropization, green line) 



but the net effect has already been described in [SI [TO]- The main features 
are shortly summarized: The thermo effect is described in the present case on 
a kinetic level by taking into account a non equilibrium distribution function 
for the plasma background, denoted by /b, in the derivation of appropriate 
Trubnikov/Rosenbluth potentials. In particular /f, has been constructed from 
a Maxwellian distribution, which is perturbed by a series of Hermitian polyno- 
mials (see |B|) . The coefficients in the perturbation series mainly depends on the 
temperature gradient and the ratio Tb/n;,. For high densities, rif, ~ 10^° /m^, 
the influence of the temperature gradient is small, because the coUisionality is 
high and the plasma is always close to local equilibrium. At low densities, e.g. 
closer to the target plates of a fusion device, even small temperature gradients 
lead to a significant temperature gradient force. Nevertheless the perturbation 
from equilibrium may not exceed the limit where /f, < 0. 

8 A reference modeling scenario for the MAST 
tokamak 

The extended EIRENE code supplemented by the trace ion module has been 
tested for a realistic modeling case for the MAST tokamak [27] . Due to the low 
aspect ratio of this tokamak {A = R/a = 1.3) neoclassical transport effects are 
enhanced, which makes MAST particularly well suited for testing the trace ion 
module. The main purpose of the simulations was to see how different impu- 
rity transport processes affect impurity profiles, which has been supplemented 
by an at least qualitative comparison of the simulation results with CCD cam- 
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Figure 4: EIRENE grid for MAST discharge #13949 
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Figure 5: CIII emission profiles from CCD camera (Mast lower divertor camera, 
by courtesy of S. Lisgo) vs. EIRENE simulation with trace ion module including 
different effects; a.) cut 1 b.) cut 2 c.) cut 3 as indicated in figure 6b; d.) radial 
temperature profile, cut 3 in figure 6b 
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era images. The present case is based on discharge #13949, during which the 
CCD system recorded CII and CIII emission. The plasma background for this 
discharge has been provided by the OSM code [H]. The quality of the OSM 
background plasma calculation is important because the emission of CII and 
CIII is strongly dependent on the local plasma conditions. The VB drift is 
acting downward for positive ions. A toroidally symmetric methane gas puff 
has been simulated, thus a 2D simulation is sufficient. The puff location was 
lower inboard at r = 0.281 m,y ~ —1.22 m. The fragmentation of the methane 
in the divertor plasma of MAST, where temperatures are about 1 — 20eV, is 
characterized by a walk through of neutral as well as ionized molecules, which 
result into the final fragmentation products C and H (and their ions) . The sin- 
gle code concept for both types of particles avoids rapid data transfer between 
codes, which would occur in the standard treatment where two kinetic codes 
are used. Moreover it is intrinsically ensured that atomic and molecular data is 
used consistently. In the present case the methane fragmentation has been de- 
scribed with the data obtained in |21| . Ionization and recombination of carbon 
atoms/ions has been described with ADAS data [TS], in particular adfl l/scd94 
and adfll/acd94- The highest charge state of carbon included in the simula- 
tions was C^"*", because the recombination rate of C^~^ is too small for being a 
significant source of C^"*" under the plasma conditions at hand. Moreover the life 
time of C'^+ in the MAST divertor is already long enough for being completely 
thermalized and a fluid description is the better and more efficient choice. The 
grid which has been used is shown in figure |4l the boundary of which has been 
specified as carbon. Physical and chemical sputtering has been turned off. 




Figure 6: a.) CII emission from CCD image, MAST discharge #13949, t = 
280ms, (by courtesy of S. Lisgo) b.) CII emission from EIRENE simulation 
with trace ion module 
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Physics effects studied with the trace ion module are friction force and ther- 
malization, thermal force, mirror force, parallel electric field force and drift 
effects. Available drifts are VU-drift, curvature drift including a quasi 3D cor- 
rection term for 2D simulations and E x B drift, which requires an appropriate 
electric field model. If no effects for ions are activated particles simply follow 
magnetic field lines. 

Radial and field aligned emission profiles with different effects active are 
shown in figure [5] and have to be understood as indicated in the 2D emission 
pattern shown in figure [7] (two radial profiles 1,2 and one field aligned 3). Five 
different cases are shown: 1.) parallel motion along field lines and friction force 
+ thermalization (CC), 2.) mirror force added (M), 3.) VB drift added (GB), 
4.) curvature drift added including quasi 3D correction (CD), 5.) E x B (FT) 
added. Profiles obtained from CCD camera data are denoted by CCD in the 
figures. The electric field has been obtained from the model described in Rcf. 

the calculation of which only involves background plasma quantities. Be- 
cause temperature gradients in the region of the radiated CIII emission (inside 
the closed flux surface region) are small, the thermal force effect is negligible in 
this case and such simulations are not shown. Temperature profiles have been 
normalized to the background plasma temperature T^, thus complete thermal- 
ization is reached at a value of 1. 

The general behavior of these emission profiles is a maximum near the sepa- 
ratrix and a decrease in opposite radial direction further in the plasma. Consid- 
ering friction force and thermalization of C'^'^ the profile is relatively flat along 
the field lines and falls off radially inward, following almost the inverse of the 
temperature profile (which mainly governs the lifetime and hence the density of 
the particles). The mirror force does not affect the radial profiles but particles 
are shifted further upstream to the inboard side along the field lines. The VB 
drift is acting in downward direction in this case and yields a decrease in den- 
sity because particles are pushed to the bottom across the separatrix or in other 
words the probability for particles crossing the separatrix is lower. The major 
part of the curvature drift is due to toroidal curvature (for which a correction 
term has been introduced, see equ. [T^ and is acting in direction of the V-B-drift. 
Thus it has a similar effect and further decreases the particle density. The last 
effect which has been added is an electric field model which slightly alters the 
density profiles to lower values. It has been observed that outside the separa- 
trix the density profiles are practically zero, where it is noted that densities and 
temperatures are considered only in regions where standard deviation is below 
10%. Whereas the emission profile is significantly influenced by the drifts the 
effect of the drifts on the temperature profile is small, see figure [5j In all cases 
C^"*" is nearly thermalized near the separatrix and the degree of thermalization 
decreases further into the plasma to 80%. Close to the separatrix the lifetime 
of is about 200 /is, which is just long enough to fully thermalize. Further 
inside the plasma the lifetime is decreasing and full thermalization cannot be 
reached anymore. The observation from the simulation point of view is that the 
temperature behavior is exclusively governed by the Coulomb collision operator 
and the drift effects mainly act on the density profiles. 
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Finally the simulation results have been compared at least qualitatively 
with uncalibrated CCD camera images for carbon emission at two wave length 
(465nm, CIII and 514nm, CII) from MAST discharge 13949, which have been 
kindly provided by S. Lisgo. It is noted that raw images of the CCD system 
(resolution ~ 5 mm) have been converted into 2D poloidal emissivity profiles 
via tomographic inversion and emission from intrinsic carbon due to sputtering 
has been subtracted. Simulated carbon density profiles have been converted 
into emission profiles using photon emissivity coefficients (PEC) according to 

^ = n,nePEC{X,n„Te) (16) 

where is the particle density of the particular species, rig and Tg are the elec- 
tron density and temperature and A is the specific wave length. The particular 
PEC's have been taken from the ADAS database {pecQQ^C-Vsu^cl.dat and 
pec96^c-vsu^c2.dat^ metastable states are not resolved). Results arc shown 
for CII emission in figure [6] and for CIII in[7l The general behavior of the simu- 




Figure 7: a.) CIII emission from CCD image, MAST discharge #13949, t = 
280ms, (by courtesy of S. Lisgo) b.) CIII emission from EIRENE simulation 
with trace ion module 

lated CII profile corresponds to the experimental one, but especially the density 
extends much further into the plasma as experimental seen. The CII emission 
has been calculated from the simulated C"*" density profile, which is governed 
mainly by the source of C"*". In this case the source of C"*" is determined by 
the break up process of methane, where finally the dissociation of CH is the 
only channel to produce C and subsequently C^. In the present simulation 
the Langer database has been used to describe the break up of the methane, 
which might allow to penetrate CH or C too far into the plasma before being 
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ionized to C"*". The influence of the present break up model has not been fully 
understood yet and will be further investigated. 

In contrast the lifetime of is much longer than for C+ and the density 
profile is more affected by the transport rather than the source profile. The 
radiation close to the X point is overestimated by the simulation, but the ra- 
diation pattern itself is close to the CCD image. This is explained by the low 
aspect ratio of MAST where the strong magnetic field curvature is significantly 
enhancing the drift effects, which govern the deposition of the particles. 

The comparison of the CCD images and the trace ion transport simulations 
with the extended EIRENE code show at least that the overall model is consis- 
tent with the experimental observations. 

Finally a comment on the computational demands of the performed simula- 
tions is given. It has been observed that approximately 10^ particle histories are 
required for the present MAST case to obtain a reasonably small standard devi- 
ation in regions of interest (< 10%). Might be different for other cases. Tracing 
10^ particles on the grid shown in figure 21 which consists of 13536 triangles, 
took approximately 10 minutes on an Intel Core2 Quad at 2,83 Ghz. The max- 
imum allowed time step for the tracing ions was 10~^s and the maximum life 
time of a particle has been limited to 0.1s. On the Opteron cluster LEO I of 
Innsbruck University an EIRENE run where exactly 10^ particle histories have 
been calculated on 96 cpu's took less than one minute. Leo I is a mixture of 
Opteron Barcelona at 2,3 Ghz and Opteron Shanghai at 2,5 Ghz. 

9 Conclusions 

A linear kinetic Monte Carlo particle transport algorithm for neutral (force- 
free) particles and linear Boltzmann collision integrals has been supplemented by 
numerical orbit integration modules as well as diffusive (Fokker-Planck) velocity 
space collision integrals. The main motivation for this development was to 
combine neutral atom transport with drift (guiding center) kinetic impurity ion 
transport in magnetized fusion edge plasmas into a single algorithm. 

The unified procedure allows to simulate both neutral and charged particles 
within a single source code and hence to avoid any numerical dissipation which 
would otherwise be caused by frequent transfer of information between neutral 
particle to charged particle modules. In particular the extended code seems to 
be well suited to deal with the rich chemistry of e.g. hydrocarbons in fusion edge 
plasmas, when a fragmentation of a single molecule proceeds via a sequence of 
neutral and charged states with quite distinct transport characteristics. Also in 
other circumstances, often in those relevant for cool dense fusion divertor plas- 
mas, such as transport of weakly ionized trace impurities near target surfaces, 
a strictly consistent transport description of neutral and charged states of an 
atom can be important. 

Clearly the optimal coordinate systems for solving transport problems for 
magnetized ions and for neutrals are at odds with each other, with the former 
necessarily being related to (inhomogencous) magnetic field properties, and the 
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latter being simply Cartesian coordinates (in which also solid surfaces and hence 
boundary conditions for the ion transport are naturally formulated). 

Wc have shown that a computationally affordable orbit integration scheme, 
which can be readily implemented into Cartesian grid based Monte Carlo solvers, 
can be constructed and is sufficiently accurate for a number of relevant appli- 
cations including those mentioned above. The boundaries of the validity range, 
for acceptable time steps (and hence acceptable CPU costs of the entire pack- 
age) are specified and found to accommodate in particular typical fusion edge 
plasma and divertor plasma conditions. 

The linear guiding center Fokkcr-Planck collision integral is integrated into 
the Monte Carlo procedure by utilizing Rosenbluth potentials, evaluated "on 
the fly" on non-Maxwellian background velocity distributions to accommodate 
kinetic thermal force effects. We have shown that this is economically possible 
and that the unavoidable singularity in these potentials at the zero of perpen- 
dicular velocity can be properly dealt with by an implicit procedure of Monte 
Carlo sampling of the post-collision ion velocities. 

A series of test cases based on semi-analytic solutions to strongly simplified 
problems has been developed to verify both the implementation of the new 
collision kernel and the explicit orbit integrator. 

As an illustrative sample of the essential capabilities we apply the extended 
code to physical and configurational model parameters chosen to simulate a 
methane gas puff experiment carried out in the compact spherical tokamak 
MAST (UK). In this magnetic configuration and plasma density range kinetic 
ion drift orbit effects can be expected to be particularly pronounced. 

Experimental results on weakly ionized carbon ion emission patterns are 
discussed with respect to relative importance of some of the new features in the 
extended code. 
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A Backward Orbit Integrator 

The algorithm for advancing particle positions in the EIRENE code is given by 

Yt+At = y* + y ■ At, 

where y is the guiding center position in case of ionized particles and y is 
an effective guiding center velocity. Computationally affordable methods for 
evaluating y, which fit to the Monte Carlo framework of EIRENE, arc Adams 
Bashforth backward methods. The backward formula for an effective guiding 
center velocity for different numbers of backward steps j is given by 

4 

8=0 

where contains coefficients corresponding to the number of backward steps 
required. If j = 1 the first order Euler method is restored, where the effective 
guiding center velocity is calculated solely from the current plasma and magnetic 
field parameters at the current location as is described in detail in chapter [5j It 
is noted that have to be integrated in the same manner as y. For a 4 step 
backward method the following coefficients have been calculated: 

" 1 

3 _i 

2 2 " 

12 12 12 
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24 24 24 
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. 720 360 60 

Guiding center orbits are started with the Euler method and take each follow- 
ing time step on more precedent guiding center velocity, until four values are 
available. These backward methods require no additional computational effort, 
just one more array has to defined which stores the values of y and v\\^± for j 
successive steps. 

B Coulomb Collision Operator for Non Max- 
wellian Plasmas 

In the test particle approach the treatment of collisions with the background 
plasma (denoted by index b), given by fluid quantities, requires the construction 
of a kinetic distribution function /{,. One approach, taking into account both 
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the friction force as well as the thermo effect, is to describe the local plasma 
equilibrium by a Maxwellian like reference state, perturbed by a distortion func- 
tion accounting for inhomogcneities. The distortion function can be a series of 
Laguerre Polynomials, e.g. |5], or Hcrmite Polynomials, e.g. [HI HH]- In the 
present work the 13 moment approximation is considered, where the deviation 
from equilibrium r]{vb, t) is given by a truncated series of Hermitian polynomials. 
This approach has been developed in [20] and applied to plasma transport codes 
in [HI dni ■ The calculation of a kinetic distribution function with this method 
requires in addition to the five fluid moments number density, nb(x), temper- 
ature, Tb(x) and flow velocity, Ub(x), three components of the heat flux q(x) 
and five independent components of the pressure tensor tt, hence 13 moments 
are used. The space coordinate is denoted by x. For deriving appropriate drift 
and diffusion coefficients in the standard Trubnikov/Rosenbluth formalism it is 
convent to use the inverse thermal velocity, vth = v'ti&/22X, as a normalization 
parameter, 

1 / mi, , , 

«= — = \/^' (17) 



and furthermore to introduce the following normalized velocities 

c = a(vh - Ub(x)), X = a(v - Ub(x)), y = c - x, (18) 

where vj, is the single particle velocity of background particles and v denotes 
the test particle velocity. Then the perturbed Maxwellian distribution function 
is given by 

/fc(x,c) = n,,(x)-^cxp(-c-) (l + ?7(x,c)). (19) 

The distortion function in the 13 moment approximation reads (summation over 
i and j) 

r7(x,c) = /.f )(x)iJ,p)(c) + h^\^)Hg\c), (20) 

with the definition of the Hermitian polynomials in terms of normalized veloci- 
ties, 

iyf)(c) = ;^c«(2c^-5) 

H^J^ic) = V2(^c,c,-icX^. (21) 

The corresponding coefficients are related to the heat flux q(x) and the pressure 
tensor tt by 
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By substitution of c = y + x the perturbed Maxwellian can be written as 

(x)q 
7r3/2 



+E'^?^.?'(y+^)| • (23) 



and the Trubnikov/Roscnbluth potentials can be written as convolution inte- 
grals of the form 

0(x,x) = / ^-.h{y + x) 

Jr^ y 



f '4'-My + x). (24) 



The potential functions can then be written as a sum of integrals depending on 
the normalized test particle velocity Xi, 

</)(x,x) = nfc(x)(/i(x)+ 

i i,j 

E^^'^^'W+E^^i'^^w) • (25) 

The integrals /^^^ and the details of the evaluation as well as the corresponding 
results are given in[Cj By introducing a generalized heat flux Q and a generalized 
pressure 11 tensor with 

and skipping ri(x) the Trubnikov/Roscnbluth potentials can be written as 

Hx) = 0o + EQ»-'^i 



+ En.(^.-3^) 

Hx) = v^o+E^'-'/'i 



En.,(5.,-3^)V'2 (27) 
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where it is noted that the potential functions </>o,i,2 and V'0.1,2, which are de- 
fined in [33l depend only on the absolute value of x = j | x 1 1 and are ab initio 
independent of the gyro angle. 

C Potential Functions of the Coulomb Collision 
Operator 

For obtaining the Trubnikov/Rosenbluth potentials for a perturbed Maxwellian 
velocity distribution function the following integrals have to be calculated 



71-^/^ Jr3 y 
^'W = 4i/ d'y^e-^y+^^'Hl'\y + x) 

lUx) -472 f d'y^ e-^y^^)'H^\y + x) (28) 

with the Hcrmitian polynomials defined as in 1211 The integration can be per- 
formed by switching to a spherical coordinate system (e^, e^, e^) with parallel 
to the vector x, then x and y have the form 

y sin 6 cos ip \ 

(29) 

and the squared absolute value of the sum of these two vectors reads 

(y + xf = x^ + y^ + 2X2/ cos 9. (30) 
In such a coordinate system integrals of the following form arise 

Hx) = ^ / / y^smOdyded^... (31) 
' Jo Jo Jo 

the integration of which is straight forward. Since the solutions of the integrals 
consist of combination of the error function and it's derivative we introduce 

$==er/(x) <S>' = ^e-^\ (32) 
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and define the following potential functions in the same way as in [5] 
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(33) 



Using these definitions one can write the solutions of the integrals in a convenient 
way 
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(34) 



These results have been calculated for a coordinate system with parallel to 
X one has to transform it back into the coordinate system of the laboratory, 
which can be done by a rotation defined by two matrices 



M,((p) 



COS if — sin (y9 
sin (y3 COS Lp 


COS 6 sin I 

— sin cos ( 



(35) 



A vector in the coordinate system (er,ee,e^) can be transformed to the labo- 
ratory system (e^;, ej,, e^) denoted by L via 



r3.4 



(e„ey,e,)M^(^)Mj,(0) 




(36) 
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and a matrix 

ll'^ = MI^'^^M^ M = R^{(p)Ry{9) (37) 

Using the following relations of the angles and the coordinates in the laboratory 
system 

— = sin 9 cos ip — = sin 9 sin p — — cos 9 (38) 
XXX 
the final results of the integrals can be written as 
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These results coincide with the ones obtained by Reiser in and |10) . 



D Derivation of the Drift and Diffusion Coeffi- 
cients 

The explicit form of the drift and diffusion coefficients |3] and the derivatives of 
the potential functions [27] is given. In contrast to [5] , where derivatives have 
been calculated for the parallel direction only, also the perpendicular direction 
is considered in the present case, thus allowing to extend the thermo effect in 
perpendicular direction. This extension is planned in a future project. The 
following helpful quantities have been used, 

■X 



dvk 



Xk. 
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r^ — \ ft 

OXl X X 



dXi 

XiXj 

X^ 



Xi_ 

X 



Then it follows for the first derivative with respect to Vk 
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X^Xk 
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(39) 
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where 



X 



The second derivative reads 

a 



dvidvk 



XkXl III 



XI 
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XiXk ,11 
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3diMikj(l)2 



3 ( —Mikj 
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X 



where the following tensors have been introduced as 



(40) 



(41) 



(42) 



XtXk 
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Xi 



X] 



diMikj = NuNkj + —diNkj + ^diN,k + N,kNij. 
In addition to the derivatives of the potential functions the scalar product 



Ci = Qx 

and the sum over all elements of the matrix 



C2 = ^ n 



,XtXj 

' 2~ 

X 



(43) 



(44) 



have to be evaluated for obtaining the explicit form of the drift and diffusion 
coefficients. The resulting scalar quantities CI and C2 are not independent 
of the gyro angle ab initio. Hence they do not fit into the drift kinetic model 
which is the aim of the present work. Neglecting perpendicular heat fluxes 
(same approach as in [5] ) the scalar product 02] is gyro angle free and can be 
evaluated by 

1 brrib 1 



CI 



VsV Jib 



where the parallel heat flux is defined as 

-,5/2 



15 T°^' le^^eg , 



(45) 



(46) 
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The value of the transport coefficient K|[ is 1/0.5657 according to [50]. With 
this approach at least the thermo effect in parallel direction can be studied. 
According to [30] this viscosity tensor is invariant under rotations around b and 
it is assumed that the scalar quantity C2 is ab initio independent of the gyro 
angle for magnetized plasmas. The reduced form of the viscosity tensor is given 

by 

12 TT^e^ 1/2 5/2 II 

"'^^"71^^"^^ ^^ ^"'^^^ ^^^^ 

where 

II 1 
" = 2 



'^Vxx 

Vyy + J/^^ 

Uuu + Vz 



Here the x direction coincides with the direction of the magnetic field vector. 
The explicit form of the components vu is 

vu = 2— Vvb 

0Xj O 



with the factor m = 1.0/(1.2 + 0.8485Z^ 



E Diffusion Coefficients for Langevin Equations 

The Langevin approach for solving Fokkcr-Planck equations requires a decom- 
position of the diffusion tensor according to D = BB^ . The choice of B is 
not unique which reflects somehow the fact that there is more information in a 
random walk trajectory than needed to solve the Fokker-Planck equation [T9] . 
In this chapter a reasonable method is given which treats the case at hand 
where the diffusion tensor is a positive definite symmetric 2x2 matrix. It can 
be decomposed as 

D = TAT^, (48) 

where T is an orthogonal matrix. Then the diagonal matrix A contains the 
eigenvalues of D. Introducing another matrix S with A = SS'^ it follows 

D = TKT'^ = TSS^T^ = TS{TSf = BB'^ . (49) 

Thus a choice for the matrix B is given by _B = TS, where S contains the square 
roots of the eigenvalues of D, which explicitly calculated read 

A± = \{Dii + D22) ± \^{Dii - D22Y + 47^12. (50) 
The matrix T can be calculated as 

T 1 ( --^12 D22-\- \ 
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where the norm of the eigenvectors is 

N = + {Dn - A+)2 = VDi2 + (I?22-A_)2. 

The expUcit form of the matrix B is then given by 



This method has been implemented into the EIRENE code for determination 
of the elements of B at each time step. 

F Semi Implicit Method for Coulomb Collisions 
of Low Energetic Particles 

Evaluating drift and diffusion coefficients [3] for purely Maxwellian distribution 
functions of the plasma background in the cylindrical coordinate system at hand 
yields a singularity in the perpendicular drift coefficient. Thus the present Monte 
Carlo procedure has to be supplemented by an implicit method for properly 
advancing discrete perpendicular velocity increments in time. The explicit form 
of the perpendicular drift coefficient written in terms of the potential function 
(j) only reads 

1 1 \ 11 



A2^11 



X V 4x2/ 2x2 



(53) 



where (j) = (j)o/a and Xi = without loss of generality, thus x = vxT+xf = X2- 
Its singular behavior near is obvious from the Taylor series expansion, 

A2{x^0)^ ^- + 0{xf. (54) 
'iV'^ X 

The normalized perpendicular velocity is formally advanced in one time accord- 
ing to 

rt+At 

^.ost^col ^ ^.re^col ^ A,{x2{t'), Xl)dt' , (55) 

which has been solved imphcitly with the first order Euler method given by 

post— col pre— col , \ j. a f post—col\ \ /rn\ 

X2 =X2 +AtA2(x2 )>Xi)- (56) 

For performing the iterations the derivative of A2 with respect to V2 is required, 
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where the coefficients a and b are defined as follows 









)-\ 







2 VX2 X^ 



The second order Newton Raphson method is sufficient to solve the actual sys- 
tem which is given by 

£/ post — col\ post — col pre— col \ _t a / post—col\ 

/(X2 ) = X2 -X2 -AtAaix^ ) 

f^^post-coi^ = 1 - At A^(xr*"™^ (59) 

and the well known iteration formulae: 

Xn+1 = Xn + . ^ , (60) 

The implicitly obtained result for ■^p"'**"^"' jg used to calculate an effective 
^2 — (X2°''*^^°' ~ X2^^~^°'^)/ which then overwrites A2 from the explicit 
scheme. 

In practice it is sufhcient to activate the presented implicit method only in 
case of small values X2 ^ 1, because the extremely stiff behavior of A2 near 
changes into a slightly increasing function for more positive values of X2- A 
threshold for X2 can be specified in the input of the code. This parameter has 
been set to 0.01 (dimensionless) for performing trace carbon simulations in the 
MAST divertor region. Actually this number depends strongly on the mass 
difference of the collision partners. If the test particle is much heavier as the 
field particle this number should be increased and vice versa. 

Moreover /(Xi°'**~'^°') is a monotonically increasing function without minima 
and maxima and less than 10 iterations are required in almost any case to reach 
a sufRcient accuracy. The iteration is currently stopped if the error of the 
normalized quantity X2, which is of the order unity, drops below 10~^. The 
additional computational costs are raised by 10% once the implicit method is 
activated, but only in this case Coulomb collisions are properly treated. 
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